Cloud Accretion Simulation: The “Warm Rain” Process#
This notebook simulates the growth of cloud droplets into raindrops through the process of accretion (also known as collision-coalescence). This mechanism is the primary driver of rain formation in “warm clouds” (clouds where temperatures are above freezing), as detailed in the course notes.
The Physics of Accretion#
In a cloud, droplets of different sizes exist simultaneously. Because larger droplets are heavier, they fall faster than smaller droplets. This difference in terminal velocity allows larger drops to overtake and collide with smaller ones in their path.
Before we start to introduce the code, let’s quickly load the modules we will need:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from IPython.display import HTML
from dataclasses import dataclass
from scipy.spatial.distance import cdist
import math
### from mpl_toolkits.mplot3d import Axes3D
# Enable interactive plots for animation in Jupyter
%matplotlib notebook
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 6
4 from IPython.display import HTML
5 from dataclasses import dataclass
----> 6 from scipy.spatial.distance import cdist
7 import math
9 ### from mpl_toolkits.mplot3d import Axes3D
10
11 # Enable interactive plots for animation in Jupyter
ModuleNotFoundError: No module named 'scipy'
1. Terminal Velocity#
A droplet’s terminal fall speed (\(w\)) depends on its radius (\(r\)).
We can simplify the relationship by allocating the behaviour to one of three distinct regimes:
Stokes’ Law (\(r < 30\mu m\)): \(w \propto r^2\). Viscous forces dominate.
Intermediate Regime (\(30\mu m < r < 1mm\)): \(w \propto r\).
Newtonian Regime (\(r > 1mm\)): \(w \propto \sqrt{r}\). Aerodynamic drag dominates.
Let’s define these relationships in the following function…
The constants X1, X2, and X3 represent the coefficients for the different fall-speed regimes.
# Physical constants for terminal velocity calculations
X1, X2, X3 = 1.2e8, 8e3, 250.
def terminal_w(r):
"""Calculation of droplet terminal fallspeed based on radius."""
# Stokes regime (small drops)
w = np.where(r < 30.e-6, X1 * r**2, X2 * r)
# Aerodynamic drag regime (large drops)
w = np.where(r > 1.e-3, X3 * np.sqrt(r), w)
return w
Now we need a function for converting radius to drop mass and back again.
This is straightforward, given the density of liquid water \(\rho_l\), the relationship is simply mass=density * volume or
\( M_{drop}=\frac{4}{3} \rho_l \pi r^3\)
rho_l = 1000. # density of liquid water (kg/m^3)
def r2mass(r):
"""Converts radius to mass."""
return (4.0 * rho_l * np.pi * r**3) / 3.0
def mass2r(m):
"""Converts mass to radius."""
r = ((3.0 * m) / (4.0 * rho_l * np.pi))**(1./3.)
return r
2. Collision and Coalescence#
When a large drop overtakes a smaller one, a collision may occur. The probability of these drops sticking together is defined by the Collision Efficiency (\(E\)).
If coalescence occurs:
The masses are combined: \(M_{new} = M_{large} + M_{small}\)
The new radius is calculated from the new mass.
The terminal velocity is updated (the drop accelerates).
The smaller drop is removed from the simulation (“dead”).
The Cloud Class#
This class manages the state of the cloud. It stores the position \((x, y, z)\) and radius \(r\) of every droplet.
The collision method is the computational core:
It calculates the pairwise distance between drops in the X-Y plane.
It checks if drops have “swapped” Z-positions within a timestep (meaning the faster one overtook the slower one).
If they overlap and overtake, a random check against efficiency \(E\) determines if they merge.
Rainfall Calculation#
In the final method of this class, we calculate the rainfall. When a drop passes below the vertical altitude \(z=0\), it is counted as rainfall and removed from the cloud.
class Cloud:
"""Class of a cloud with arrays of x, y, z, r, and w."""
def __init__(self, xpos, ypos, zpos, radius, collision_E):
self.x = xpos
self.y = ypos
self.z = zpos
self.r = radius
self.dead = np.zeros(len(radius), dtype=bool)
self.rain = 0.0
self.collision_E = collision_E
# Initialize vertical velocity based on radius
self.w = terminal_w(radius)
self.znew = np.copy(self.z)
def collision(self):
ndrop = len(self.r)
# Calculate pairwise sum of radii (collision threshold)
# We use a trick to avoid a double loop for speed, though it is memory intensive for large N
rsq = np.tile(self.r, (ndrop, 1))
rsum = rsq + np.transpose(rsq)
# 1. Check Horizontal Proximity
coords = np.column_stack((self.x, self.y))
xydist = cdist(coords, coords)
# 2. Check Vertical Overtake
# Did they swap z-locations during this timestep?
zd1 = np.subtract.outer(self.z, self.z)
zd2 = np.subtract.outer(self.znew, self.znew)
# Overtake condition: Relative vertical distance changed sign
# Collision condition: XY dist < radius sum AND overtake happened
ind1, ind2 = np.where((xydist < rsum) & (zd1 * zd2 < 0.0))
# Filter unique pairs and ignore self-collision
unique = (ind1 < ind2)
ind1 = ind1[unique]
ind2 = ind2[unique]
# Process collisions
for i1, i2 in zip(ind1, ind2):
if self.dead[i1] or self.dead[i2]:
continue
# Stochastic collision efficiency check
if np.random.uniform(low=0, high=1) > self.collision_E:
continue
# CONSERVATION OF MASS: Merge 2 into 1
mass1 = r2mass(self.r[i1])
mass2 = r2mass(self.r[i2])
# Update Drop 1 (The Collector)
self.r[i1] = mass2r(mass1 + mass2)
self.w[i1] = terminal_w(self.r[i1])
self.znew[i1] = min(self.znew[i1], self.znew[i2]) # Take the lower position
# Update Drop 2 ( The Collected / Dead)
self.dead[i2] = True
self.r[i2] = self.w[i2] = 0.0
self.znew[i2] = self.z[i2] = -1.e6 # Move well below ground
def raincalc(self, xysize):
"""Check for drops reaching the surface (z < 0)."""
rain_idx = np.argwhere(self.z < 0.0)
# Calculate mass hitting the ground
mass_vec = r2mass(self.r[rain_idx])
rain_amount = np.sum(mass_vec) / xysize
self.rain += rain_amount
# Remove rain drops from simulation
self.dead[rain_idx] = True
self.r[rain_idx] = 0.0
self.w[rain_idx] = 0.0
Part 3: Simulation Parameters & Initialization#
Here we set the size of the domain, the time steps, and the initial population of droplets in a python configuration class, a nice way to handle simulation options as you can also specify the option type.
Experiment Parameters:
bimodal: IfTrue, we insert a few large drops manually to kickstart accretion. IfFalse, we use a Gamma distribution (typical for real clouds).ccn: Cloud Condensation Nuclei concentration. Higher values = more, smaller drops (polluted air).collision_E: Efficiency. 1.0 means every collision results in a merger.
# ===========================
# CONFIGURATION CLASS
# ============================
@dataclass
class SimOptions:
"""Holds all parameters for the simulation."""
# Domain & Time
xsize: float = 0.0001
ysize: float = 0.0001
zsize: float = 2000.0
dt: float = 60.0
lentime: int = 180
# Physics
zcloud1: float = 1500.0
zcloud2: float = 2000.0
collision_E: float = 0.3
ccn: float = 200.0
rsmall: float = 10.e-6
rsigma: float = 5.e-6
# Distributions
bimodal: bool = False
rlarge: float = 25.e-6
ratio_large: float = 0.05
# Visualization Optimization
slice_step: int = 2 # Draw only every Nth drop
steps_per_frame: int = 2 # Run physics N times per video frame
@property
def xysize(self):
return self.xsize * self.ysize
# ==========================================
# INITIALIZATION ROUTINE
# ==========================================
def initialize_simulation(opts: SimOptions):
"""Sets up the drops and creates the Cloud object."""
ccn_m3 = opts.ccn * 1.e6
ndrop = int(ccn_m3 * opts.xysize * (opts.zcloud2 - opts.zcloud1))
print(f"Initializing Cloud: {ndrop} drops (CCN={opts.ccn})")
if opts.bimodal:
nlarge = int(np.ceil(opts.ratio_large * ndrop))
dropr = np.full(ndrop, fill_value=opts.rsmall)
dropr[-nlarge:] = opts.rlarge
print(f"- Mode: Bimodal ({nlarge} large drops added)")
else:
gshape = (opts.rsmall / opts.rsigma)**2
gscale = opts.rsmall / gshape
dropr = np.random.gamma(gshape, scale=gscale, size=ndrop)
print(f"- Mode: Gamma Distribution")
# Stats for plotting
stats = {
'rmean': np.mean(dropr) * 1e6,
'rmax': np.max(dropr) * 1e6,
'rainavail': np.sum(r2mass(dropr)) / opts.xysize,
'ndrop': ndrop
}
# Positions
dropx = np.random.uniform(0, opts.xsize, ndrop)
dropy = np.random.uniform(0, opts.ysize, ndrop)
dropz = np.random.uniform(opts.zcloud1, opts.zcloud2, ndrop)
cloud = Cloud(dropx, dropy, dropz, dropr, opts.collision_E)
return cloud, stats
def run_visual_simulation(cloud, opts: SimOptions, stats):
"""Runs the loop and generates the animation."""
print("Starting Physics & Animation Loop...")
%matplotlib inline
plt.rcParams['animation.embed_limit'] = 50.0
fig = plt.figure(figsize=(8, 8))
# Setup Axes
ax0 = fig.add_subplot(211, projection='3d')
ax0.set_xlim(0, opts.xsize * 1000.)
ax0.set_ylim(0, opts.ysize * 1000.)
ax0.set_zlim(0, opts.zsize)
ax0.set_title(f'Cloud Sim (CCN={opts.ccn})')
ax0.set_xlabel("x (mm)"); ax0.set_ylabel("y (mm)"); ax0.set_zlabel("z (m)")
ax1 = fig.add_subplot(212)
ax1.set_xlim(0, opts.lentime)
ax1.set_ylim(0, math.ceil(stats['rainavail']))
ax1.set_xlabel("Time (mins)"); ax1.set_ylabel("Rainfall (mm)")
# Init Plots
line, = ax1.plot([], [], 'r.', ms=5)
ax1.plot([0, opts.lentime], [stats['rainavail'], stats['rainavail']], 'k--', label="Max Rain")
ax1.legend(loc='upper right')
info = f"Mean R: {stats['rmean']:.1f}um\nMax R: {stats['rmax']:.1f}um\nE: {opts.collision_E}"
ax1.text(2, stats['rainavail']/2, info, ha='left')
sc3d = ax0.scatter(cloud.x[::opts.slice_step]*1000.,
cloud.y[::opts.slice_step]*1000.,
cloud.z[::opts.slice_step],
c='blue', s=cloud.r[::opts.slice_step]*1.e6)
# State containers
history = {'times': [0], 'rains': [0.0]}
sim_state = {'time': 0.0}
def animate(frame):
# Physics Sub-stepping
for _ in range(opts.steps_per_frame):
cloud.znew = np.copy(cloud.z) - opts.dt * cloud.w
cloud.collision()
cloud.z = np.copy(cloud.znew)
cloud.raincalc(opts.xysize)
sim_state['time'] += opts.dt
# Update Graphics
history['rains'].append(cloud.rain)
history['times'].append(sim_state['time'] / 60.)
sc3d._offsets3d = (cloud.x[::opts.slice_step]*1000.,
cloud.y[::opts.slice_step]*1000.,
cloud.z[::opts.slice_step])
sc3d._sizes = (cloud.r[::opts.slice_step]*1.e6)
live = int(stats['ndrop'] - np.sum(cloud.dead))
ax0.set_title(f"Time={int(sim_state['time']/60)}m, Live Drops={live}")
line.set_data(history['times'], history['rains'])
return sc3d, line
frames = int(opts.lentime * 60 / (opts.dt * opts.steps_per_frame))
ani = animation.FuncAnimation(fig, animate, frames=frames, interval=100, blit=False)
plt.close() # Prevent static plot
# KEY FIX: Use display() instead of returning the HTML object
display(HTML(ani.to_jshtml()))
# return HTML(ani.to_jshtml())
def run_cloud_experiment(**kwargs):
"""
Main entry point.
Usage: run_cloud_experiment(ccn=500, bimodal=True, zcloud1=1000)
"""
# 1. Setup Options (Override defaults with user arguments)
opts = SimOptions(**kwargs)
# --- PRINT CONFIGURATION REPORT ---
print("="*40)
print(f"{'CLOUD SIMULATION CONFIGURATION':^40}")
print("="*40)
for key, value in opts.__dict__.items():
print(f"{key:<20} : {value}")
print("-" * 40)
# 2. Initialize
cloud, stats = initialize_simulation(opts)
# 3. Run & Animate
run_visual_simulation(cloud, opts, stats)
return
Part 4: Run Simulation & Animate#
The code below creates a 3D animation.
Top Panel: The 3D view of the cloud. The size of the blue dots represents the droplet radius.
Bottom Panel: The accumulated rainfall over time.
Watch how the “rain” line stays flat initially (while drops are small and falling slowly), then spikes upwards as larger drops form via accretion and fall out quickly.
run_cloud_experiment()
========================================
CLOUD SIMULATION CONFIGURATION
========================================
xsize : 0.0001
ysize : 0.0001
zsize : 2000.0
dt : 60.0
lentime : 180
zcloud1 : 1500.0
zcloud2 : 2000.0
collision_E : 0.3
ccn : 200.0
rsmall : 1e-05
rsigma : 5e-06
bimodal : False
rlarge : 2.5e-05
ratio_large : 0.05
slice_step : 2
steps_per_frame : 2
----------------------------------------
Initializing Cloud: 1000 drops (CCN=200.0)
- Mode: Gamma Distribution
Starting Physics & Animation Loop...
It is now super straightforward to run sensitivity tests…
# An example sensitivity test....
run_cloud_experiment(ccn=600, bimodal=True, collision_E=0.5)
========================================
CLOUD SIMULATION CONFIGURATION
========================================
xsize : 0.0001
ysize : 0.0001
zsize : 2000.0
dt : 60.0
lentime : 180
zcloud1 : 1500.0
zcloud2 : 2000.0
collision_E : 0.5
ccn : 600
rsmall : 1e-05
rsigma : 5e-06
bimodal : True
rlarge : 2.5e-05
ratio_large : 0.05
slice_step : 2
steps_per_frame : 2
----------------------------------------
Initializing Cloud: 3000 drops (CCN=600)
- Mode: Bimodal (150 large drops added)
Starting Physics & Animation Loop...